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We present what we believe to be the first formal verification of a biologically 
realistic (nonlinear ODE) model of a neural circuit in a multicellular organism: 
Tap Withdrawal (TW) in C. Elegans , the common roundworm. TW is a reflex¬ 
ive behavior exhibited by C. Elegans in response to vibrating the surface on 
which it is moving; the neural circuit underlying this response is the subject 
of this investigation. Specifically, we perform reachability analysis on the TW 
circuit model of Wicks et al. (1996), which enables us to estimate key circuit 
parameters. Underlying our approach is the use of Fan and Mitra’s recently de¬ 
veloped technique for automatically computing local discrepancy (convergence 
and divergence rates) of general nonlinear systems. We show that the results we 
obtain are in agreement with the experimental results of Wicks et al. (1995). As 
opposed to the fixed parameters found in most biological models, which can only 
produce the predominant behavior, our techniques characterize ranges of param¬ 
eters that produce (and do not produce) all three observed behaviors: reversal 
of movement, acceleration, and lack of response. 

1 Introduction 

Although neurology and brain modeling/simulation is a popular field of biologi¬ 
cal study, formal verification has yet to take root. There has been cursory study 
into neurological model checking (see Section [2]), but not with the nonlinear 
ODE models used by biologists. We believe that the insights gained through 
formal verification and analysis can transform the field, as has been the case in 
the Electronic Design Automation (EDA) industry, which is now valued at over 
$4 billion annually. As EDA has allowed for increased complexity for a smaller 
time investment in hardware circuits, we believe that the same kind of benefits 
can be realized for neural circuits. 

For our initial neurological study, we have selected the round worm, Caenorhab- 
ditis Elegans, due to the simplicity of its nervous system (302 neurons, ~5,000 
synapses) and the breadth of research on the animal. The complete connectome 
of the worm is documented, and there have been a number of interesting exper¬ 
iments on its response to stimuli. 

For model-checking purposes, we were particularly interested in the tap with¬ 
drawal (TW) neural circuit. The TW circuit governs the reactionary motion of 
the animal when the petri dish in which it swims is perturbed. (A related circuit, 
touch sensitivity , controls the reaction of the worm when a stimulus is applied to 


a single point on the body.) Studies of the TW circuit have traditionally involved 
using lasers to ablate the different neurons in the circuit of multiple animals and 
measuring the results when stimuli are applied. 

A model of the TW circuit was presented by Wicks, Roehrig, and Rankin 
in EH3- This model is in the form of a system of nonlinear ODEs, as well as 
mapped polarities of the various neurons involved in the TW circuit. Addition¬ 
ally, Wicks and Rankin had a previous paper in which they measure the three 
possible reactions of the animals to TW with various neurons ablated m ; see 
also Fig. [2] The three behaviors—acceleration, reversal of movement, and no 
response—are logged with the percentage of the experimental population to dis¬ 
play that behavior. 

The |16| model has a number of circuit parameters, such as gap-junction 
conductance, capacitance, and leakage current, that crucially affect the behav¬ 
ior of the organism. Values for these parameters based on estimates, rules of 
thumb and measurements are given in m, but not the parameter ranges. A 
quick analysis of the circuit shows that variations in these parameter values may 
give rise to changes in the behavior of the model from acceleration or rever¬ 
sal to no-response. As all biological parameters vary across populations, time, 
and environments, identifying parameter ranges corresponding to behaviors is a 
fundamentally important problem. For a complete characterization of the TW 
circuit, it is therefore critical to identify the range of parameter values that give 
rise to these different types of behavior. 

Using automatically generated local discrepancy functions 03], we are able 
to perform reachability analysis on the m model. This approach combines static 
analysis with numerical simulations to allow us to iteratively compute more pre¬ 
cise over-approximations of the reachable states of the system with respect to 
a continuous range of parameter values. We used this approach, which we refer 
to as 8-refinement , to determine parameter ranges that produce all three behav¬ 
iors for the control group (no ablation) and four ablation experiments. This is a 
significant expansion of the biological results, where only static parameters are 
obtained, and only for one behavior per experimental group. The specific pa¬ 
rameters of interest are the gap-junction conductances for three sensory neurons 
in the TW circuit, as the gap junctions formed by these neurons are known to 
be the most important functional connections to the TW process. 

Our results are further organized into how many of these conductances we 
considered simultaneously, a categorization we refer to as 1-D, 2-D, and 3-D. For 
the 1-D and 2-D cases, we were able to determine parameter ranges for which 
the three TW responses are guaranteed to hold. The 3-D case is only applicable 
to the control group; here, again, we were able to produce the same guarantees. 
Moreover, with a single exception (which Wicks himself has experienced), our 
results match the trends (in terms of relative percentages) of the earlier biological 
experiments (see Fig. [2]) . 

The rest of the paper develops along the following lines. Section [3] provides 
requisite background material on the TW neural circuit, its reactionary behavior, 
and the ODE model of m • Section [4] describes our reach-tube reachability anal- 


ysis and associated property checking. Section[5]presents our extensive collection 
of model-checking/parameter-estimation results. Section [2] reviews related work. 
Section [6] offers our concluding remarks and directions for future work. 

2 Related Work 

Iyengar et al. m present a Pathway Logic (PL) model of neural circuits in the 
marine mollusk Aplysia. Specifically, the circuits they focus on are those involved 
in neural plasticity and memory formation. PL systems do not use differential 
equations, favoring qualitative symbolic models. They do not argue that they 
can replace traditional ODE systems, but rather that their qualitative insights 
can support the quantitative analysis of such systems. Neurons are expressed 
in terms of rewrite rules and data types. Using the PL formalism, they are 
able to simulate neural circuits and perform qualitative in silico experiments, 
such as simulating knock-out of an individual components or other changes to 
the network. Their simulations, unlike our reachability analysis, do not provide 
exhaustive exploration of the state space. Additionally, PL models are abstrac¬ 
tions usually made in collaboration between computer scientists and biologists. 
Our work meets the biologists on their own terms, using the pre-existing ODE 
systems developed from physiological experiments. 

Tiwari and Talcott m build a discrete symbolic model of the neural cir¬ 
cuit Central Pattern Generator (CPG) in Aplysia. The CPG governs rhythmic 
foregut motion as the mollusk feeds. Working from a physiological (non-linear 
ODE) model, they abstract to a discrete system and use the Symbolic Analy¬ 
sis Laboratory (SAL) model checker to verify various properties of this system. 
They cite the complexity of the original model and the difficulty of parameter es¬ 
timation as motivation for their abstraction. Their discrete model synchronously 
composes 10 input/output automaton (neurons), connects them with 3 types of 
links (excitatory synapse, inhibitory synapse, gap-junction), and includes an ob¬ 
server component. The input of each neuron can be positive, negative, or zero 
and the output is a boolean: a pulse is generated or not. Our approach uses 
the original biological model of the TW circuit of C. Elegans nu, and through 
reachability analysis, we obtain the parameter ranges of interest. 

We have extensive experience with model checking and reachability analysis 
in the cardiac domain, e.g. j7l9linl13( . In fact, much of our previous work has 
focused on the cardiac myocyte, a computationally similar cell to the neuron. 
This is not surprising as both belong to the class of excitable cells. The similarities 
are so numerous that we have used a variation of the Hodgkin-Huxley model of 
the squid giant axon [8] to model ion channel flow in cardiac tissue. 

3 Background 

In C. Elegans , there are three classes of neurons: sensory , inter , and motor. For 
the TW circuit, the sensory neurons are PLM , PVD , ALM, and AVM , and the 




inter-neurons are AVD, DVA, PVC, AVA, and AVB. The model we are using 
abstracts away the motor neurons as simply forward and reverse movement. 

Neurons are connected in two ways: electrically via bi-directional gap junc¬ 
tions, and chemically via uni-directional chemical synapses. Each connection has 
varying degrees of throughput, and each neuron can be excitatory or inhibitory, 
governing the polarity of transmitted signals. These polarities were experimen¬ 
tally determined in [16], and used to produce the circuit shown in Fig. [l] 



Fig. 1 . Tap Withdrawal Circuit of C. Elegans. Rectangle: Sensory Neurons; Circle: 
Inter-neurons; Dashed Undirected Edge: Gap Junction; Solid Directed Edge: Chemical 
Synapse; Edge Label: Number of Connections; Dark Gray: Excitatory Neuron; Light 
Gray: Inhibitory Neuron; White: Unknown Polarity. FWD: Forward Motor system; 
REV: Reverse Motor System. 


The TW circuit produces three distinct locomotive behaviors: acceleration, 
reversal of movement, and a lack of response. In m , Wicks et al. performed 
a series of laser ablation experiments in which they knocked out a neuron in a 
group of animals (worms), subjected them to a tapped surface, and recorded 
the magnitude and direction of the resulting behavior. In the control group 
with no neurons knocked out, 98% of subjects reacted to a tap with a reversal of 
locomotion, but there were still measured cases of acceleration and “no response” 
behavior. Fig. [2] shows the response types for each of their experiments. 

The dynamics of a neuron’s membrane potential, V, is determined by the 
sum of all input currents, written as: 

C m V = -^-(V t -V) + Yi J9ap + isyn + istim 
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Fig. 2. Effect of ablation on Tap Withdrawal reflex. The length of the bars indicate 
the fraction of the population demonstrating the particular behavior. 


where C m is the membrane capacitance , R m is the membrane resistance , V) is 
the leakage potential , I 9ap and I syn are gap-junction and the chemical synapse 
currents, respectively, and I stlm is the applied external stimulus current. The 
summations are over all neurons with which this neuron has a (gap-j unction or 
synaptic) connection. 

The current flow between neuron i and j via a gap-junction is given by: 

nr = <; p g r(Vj - vo 

where the constant g^ p is the maximum conductance of the gap junction, and 
n ij P the number of gap-junction connections between neurons i and j. The 
conductance g 9 ^ 9 is one of the key circuit parameters of this model that dra¬ 
matically affects the behavior of the animal. 

The synaptic current flowing from pre-synaptic neuron j to post-synaptic 
neuron i is described as follows: 





















where gij' l (t) is the time-varying synaptic conductance of neuron i, n s dj n is the 
number of synaptic connections from neuron j to neuron i, and Ej is the reversal 
potential of neuron j for the synaptic conductance. 

The chemical synapse is characterized by a synaptic sign, or polarity, spec¬ 
ifying if said synapse is excitatory or inhibitory. The value of Ej is assumed 
to be constant for the same synaptic sign; its value is higher if the synapse is 
excitatory rather than inhibitory. 

Synaptic conductance is dependent only upon the membrane potential of 
presynaptic neuron Vy, given by: 


9TU=9T(Vj) 

where g^ n is the steady-state post-synaptic conductance in response to a pre¬ 
synaptic membrane potential. 

The steady-state post-synaptic membrane conductance is modeled as: 
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where g^( n is the maximum post-synaptic membrane conductance for the synapse, 
VeQj is the pre-synaptic equilibrium potential , and VRange is the pre-synaptic 
voltage range over which the synapse is activated. 

Combining all of the above pieces, the mathematical model of the TW circuit 
is a system of nonlinear ODEs, with each state variable defined as the membrane 
potential ofa neuron in the circuit. Consider a circuit with N neurons. The 
dynamics of the i th neuron of the circuit is given by: 
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The equilibrium potentials (Veq) of the neurons are computed by setting the 
left-hand side of Eq. ([lj to zero. This leads to a system of linear equations, that 
can be solved as follows: 


Veq = A x b 


where matrix A is given by: 


( 5 ) 







and vector b is written as: 


N 

b i =V li +R mi '£ l E j n% n gT/2- 

j=l 

The potential of the motor neurons AVB and AVA determine the observable 
behavior of the animal. If the integral of the difference between Vava ~ Vavb is 
large, the animal will reverse movement. By extension, if the difference is a large 
negative value, the animal will accelerate, and if the difference is close to zero 
there will be no response. The equation that converts the membrane potential 
of AVB and AVA to a behavioral property, (e.g. reversal ), is given by: 

Propensity to Reverse cc J(Vava — VA V B)dt (6) 

where the integration is computed from the beginning of tap stimulation until 
either the simulation ends or the integrand changes sign. To allow initial tran¬ 
sients after the tap, the test for a change of integrand sign occurs only after a 
grace period of 100 ms. 

For the purpose of reachability analysis (Section [d]), we normalize the system 
of equations with respect to the capacitance. Combining Eqs.( [lj and ( |4j) and 
taking C mi to the right-hand side, we have: 
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Now letting g\ eak = R J Cm , gf ap = fgg, and If xt = ^ the 

system dynamics can be written as: 
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This is the 9 dimensional ODE model of the TW circuit. The key circuit param¬ 
eters are the gap conductances, gf ap , and we aim to characterize the ranges of 
these conductances that produce acceleration, reversal, and no response. 


4 Reachability Analysis of Nonlinear TW Circuit 

Reachability analysis for verifying properties for general nonlinear dynamical 
systems is a well-known hard problem. Our approach relies on a recent line of 
investigation that combines static analysis of the model with validated numerical 
simulations mm- 











4.1 Background on Reachability using Discrepancy 

Consider an n-dimensional autonomous dynamical system: 

x = f{x), (8) 

where / : R™ —> R" is a Lipschitz continuous function. A solution or a trajectory 
of the system is a function £ : R™ x R>o —>• R n such that for any initial point 
G R" and at any time t > 0, f(xo,t) satisfies the differential equation Q. 
A state x in R" is reachable from, the initial set 0 C R ra within a time interval 
[ti, £ 2 ] if there exists an initial state xq £ 0 and a time f £ \t\ , £ 2 ] such that 
x = £(xo,t). The set of all reachable states in the interval [ti,t 2 ] is denoted 
by Reach(6>, [ti,^ 2 ])- If U = 0, we write Reach(t 2 ) when set 0 is clear from 
the context. If we can compute or approximate the reach set of such a model, 
then we can check for invariant or temporal properties of the model. Specifically, 
C. Elegans TW properties such as accelerated forward movement or reversal of 
movement fall into these categories. Our core reachability algorithm m uses 
a simulation engine that gives sampled numerical simulations of Q. 

Definition 1. A (xq, t, e, T)-simulation of a is a sequence of time-stamped 
sets (Ro,t 0 ), (Ri, h )..., (R n , t n ) satisfying: 

1. Each Ri is a compact set in R” with dia(Ri) < e. 

2. The last time t n =T and for each i, 0 < ti — 1 < r, where the parameter 

t is called the sampling period. 

3. For each ti, the trajectory from Xq at ti is in Ri, i.e., £(xo,ti) G Ri, and for 

any t £ the solution £(xo,t) £ hull(Ri-\, Ri). 

The algorithm for reachability analysis uses a key property of the model called 
a discrepancy function. 

Definition 2. A uniformly continuous function j3 : R n x R™ x R>o —> R>o is a 
discrepancy function of 0 if 

1. for any pair of states x, x' £ R ra , and any time t > 0, 

U(x,t) - < /3(x,x',t), and (9) 

2. for any t, as x —► x', f) —> 0. 

If a function /? meets the two conditions for any pair of states x, x' in a compact 
set K then it is called a K -local discrepancy function. Uniform continuity means 
that Ve > 0, Vx, x' £ K, 35 such that for any time t, ||x—x'|| < 6 => j3(x, x', t) < e. 
The verification results in EEEIII required the user to provide the discrepancy 
function /3 as an additional input for the model. A Lipschitz constant of the 
dynamic function / gives an exponentially growing /3, contraction metrics m 
can give tighter bounds for incrementally stable models, and sensitivity analysis 
gives tight bounds for linear systems 2] ; but none of these give an algorithm for 
computing /3 for general nonlinear models. Therefore, finding the discrepancy 
can be a barrier in the verification of large models like the TW circuit. 






Here, we use Fan and Mitra’s recently developed approach that automat¬ 
ically computes local discrepancy along individual trajectories [8]. Using the 
simulations and discrepancy, the reachability algorithm for checking properties 
proceeds as follows: Let the U be the set of states that violate the invariant in 
question. First, a (5-cover C of the initial set 0 is computed; that is, the union 
of all the (5-balls around the points in C contain 0. This S is chosen to be large 
enough so that the cardinality of C is small. Then the algorithm iteratively and 
selectively refines C and computes more and more precise over-approximations 
of Reach(@,T) as a union U Xo gcReach T). Here, Reach(£?, 5 (:ro), T) is 
computed by first generating a (xq, t, e, T)-simulation and then bloating it by 
a factor that maximizes /3(x,x',t) over x,x' £ Bs(sq) and t £ [t*_If 
Reach (Bs{xq), T) is disjoint from U or is (partly) contained in U, then the al¬ 
gorithm decides that Bg(x 0 ) satisfies and violates U, respectively. Otherwise, a 
finer cover of B$(x o) is added to C and the iterative selective refinement contin¬ 
ues. We refer to this in this paper as (5-refinement. In [3j, it is shown that this 
algorithm is sound and relatively complete for proving bounded time invariants. 

4.2 Applying Local Discrepancy to TW Circuit 

Fan and Mitra’s algorithm (see details in 0) for automatically computing local 
discrepancy relies on the Lipschitz constant and the Jacobian of the dynamic 
function, along with simulations. The Lipschitz constant is used to construct 
a coarse, one-step over-approximation S of the reach set of the system along 
a simulation. Then the algorithm computes an upper bound on the maximum 
eigenvalue of the symmetric part of the Jacobian over S, using a theorem from 
matrix perturbation theory. This gives a piecewise exponential /?, but the expo¬ 
nents are tight as they are obtained from the maximum eigenvalue of the linear 
approximation of the system in S. This means that for models with convergent 
trajectories, the exponent of /3 over S will be negative, and the Reach (T) ap¬ 
proximation will quickly become very accurate. In the rest of this section, we 
describe key steps involved in making this approach work with the TW circuit. 

The model of the TW circuit from Section [ 3 ] can be written as V = f(V), 
where V £ M 9 . The Jacobian of the system is the matrix of partial derivatives 
with the ij th term given by: 
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For parameter-range estimation of the TW circuit, each parameter p of in¬ 
terest is added as a new variable with constant dynamics (p = 0). Computing 
the reach-set from initial values of p is then used to verify or falsify invariant 
properties for a continuous range of parameter values, and therefore a whole 
family of models, instead of analyzing just a single member of that family. Here 








the parameters of interest are the quantities 1 / g \ eak , 10/ g 9 i ap , 1 / g^ yn . Consider, 
for example, 1 /g[ eak as a parameter: 
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l/g Ua.k 




In this case the Jacobian matrices for the system with parameters will be singular 
because of the all-zero rows that come from the parameter dynamics. The zero 
eigenvalues of these singular matrices are taken into account automatically by 
the algorithm for computing local discrepancy. In this paper we focus on 10 /g yap , 
leaving the others for future work. 


4.3 Checking Properties 

Once the reach sets are computed, checking the acceleration , reversal , and no¬ 
response properties are conceptually straightforward. For instance, Equation ©> 
gives a method to check reversal movement. Instead of computing the integral 
of (Vava ~ Vavb), we use the following sufficient condition to check it: 

V t G T mt ,y x G Reach(0, [t,t)),V A vA{x) > V AV b{x)- 

Here, is a specific time interval after the stimulation time, 0 is the initial set 
with parameter ranges, and recall that Reach (0, [f, f]) is the set of states reached 
at time t from 0. We implement this check by scanning the entire reachtube 
and checking that its projection on Vavb{x) is above that of Vava {x) over all 
intervals. If this check succeeds (as in Figure |4ja)), we conclude that the range 
of parameter values produce the reversal movement. If the check fails, then the 
reversal movement is not provably satisfied (Figure [4](b) ) and in that case we 
(5-refine the initial partition. 
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(a) Rev. property satisfied with g^y M ~ 1000- 


(b) Rev. unknown with g 9 ^y M = 33.33. 


Fig. 3. Model Checking Reversal Property of Control Group, with <5 = 5e — 5, varying g^y^ M - 
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(b) Rev. property satisfied with <5 = 5e-5. 


Fig. 4. Model Checking Reversal Property of Control Group by refining <5. 


5 Experimental Results 

In this section, we apply reachability analysis to parameter rangers that produce 
three different behaviors ( reversal, acceleration, no response) in the control and 
four ablation groups. When a tap stimulus is applied, the sensory neurons ( ALM , 
AVM and PLM ) propagate that signal to the motor neurons via interneurons. 
The gap-junctions formed by the sensory neurons are the most important func¬ 
tional connections to this process pQ. Therefore we vary only the gap-junction 
conductance, g pap , of the sensory neurons and keep all other parameters con¬ 
stant, as per nsj. Our experiments characterize parameter ranges for reversal, 
acceleration and no response behaviors in all groups, except the ALM,AVM- 
group where reversal behavior is not observed. 

In section |4j we explain that we use 10/g 9ap as our parameter in the state 
vector instead of g pap . Assume pf ap = 10/g? ap , i £ {AVM, ALM, PLM}. The 
corresponding range for p pap is [0.01,1]. From the reachability analysis, we esti¬ 
mate ranges for p pap that can be converted back to gf ap . 

In the following subsections, we will present our results for parameter range 
estimation for all three behaviors of the control and ablation groups. This process 
requires three experiments per group. 

5.1 1-D Parameter Space 

Here we vary one conductance at a time for two groups: the control and the 
ALM, AVM ablation groups. 

Control: For the control case, we varied P 9 avm- We found that the reversal 
property is satisfied in sub-range [0.01, 0.214] with 5 = le — 6 and the accelera¬ 
tion property in sub-range [0.63,1] with 5 = le — 5. Recall, S is the size of the 
finest cover used by the verification algorithm. We could not verify any prop¬ 
erty for the sample points in sub-range (0.214,0.63). As shown in Table [l] the 
















parameter range producing reversal , as identified by our procedure, dominates 
the parameter range for acceleration. Our procedure also shows that no value of 
P 9 avm produces the no-response behavior for the control group. 

The time required for our procedure is dependent upon the property, the 
interval for each dimension, and the size of S. For example, the time necessary 
to complete the procedure for the reversal property is approximately one hour. 


ALM, AVM Ablation Group: In this group two sensory neurons, ALM , and 
AVM , are ablated. As such, we vary only p 9 p[ M - Acceleration is satisfied over 
the interval [0.01, 0.3] with 5 = 5e — 5 and no response behavior is satisfied over 
[0.75,1] with the same S. Despite using a very small S for refinement, we did not 
observe any reversal behavior in this entire range. Examining Table [T] we see 
that acceleration is the dominant behavior for this group. 

5.2 2-D Parameter Space 

We lead with results for the control group, then examine various ablation groups. 

Parameter Refinement in 2-D: Fig. [5] helps paint a picture of how the 6- 
refinement process discussed in Section [4] works in 2-D. We consider 4 refinement 
steps for the control group: d = 7e — 5, 6 = 6e — 5, 6 = 5.5e — 5, and S = 5e — 5. 
For S = 7e — 5, the property of interest is unknown at all points. With S = 6e — 5 
the property is considered unknown for all red areas in the figure, including red 
and blue areas. Blue areas show where S = 5.5e — 5 are satisfied, and in the blue 
and yellow area both <5 = 6e — 5 and 6 = 5.5e — 5 have a satisfied property. The 
property is satisfied for the entire range of the graph when <5 = 5e — 5. Thus, the 
refinement process stops at 5 = 5e — 5, and the entire range of the 2-D parameter 
space is characterized. 


Control Here we consider the p 9 avm an d P 9 alm conductances simultaneously. 
For this group, reversal is satisfied over the range [0.01,0.0105] with S = 2e — 5 
and acceleration is satisfied over [0.63,0.6305] with the same S. Table [l] shows 
that reversal , like in the 1-D case, dominates and no response is not generated. 


PLM Ablation Group: As the PLM neuron is ablated in this group, varying 
only p 9 avm an d P 9 alm sufficient to produce all three behaviors. Here we find re¬ 
versal satisfied over [0.01,0.0105] with S = 2e —5, acceleration over [0.67,0.6705] 
with S = 5e — 5, and no response over [0.9995,1] with <5 = 5e — 5. Table [I] shows 
that reversal dominates the other two behaviors, but all three are produced. 


ALM Ablation Group: To produce all three behaviors of this group we vary 
only p 9 avm an d P 9 plm ■ Reversal is satisfied over the interval [0.01,0.0105] with 
S = 5e — 5, acceleration over [0.67, 0.6705] with 5 = 2e — 5, and no response over 



Fig. 5. Example of 2-D Parameter Refinement. Red Regions are Unknown for both 
S = 6e — 5 and $ = 5.5e — 5, Red/Blue Regions are Unknown for 5 = 6e — 5, but 
Satisfied for <5 = 5.5e — 5, and Yellow/Blue Regions are Satisfied for both. 


Group Name 

Case 

Normalized gap-junction conductance (gj gap ) 

Rev. 

Ace. 

No Resp. 

Control 

ID 

953.21 

5.873 

0 

2D 

2267.57 

1.58e-4 

0 

3D 

970.5931 

2e-6 

1.28e-7 

PLM - 

2D 

2267.57 

1.24e-4 

2.503e-7 

ALM- 

2D 

2267.57 

1.24e-4 

2.503e-7 

ALM,DVA- 

2D 

148.72 

1.24e-4 

2.503e-7 

ALM,AVM- 

ID 

0 

966.67 

3.33 


Table 1. Regions in the parameter space in which the properties are proven satisfied. 
A 1-D case shows interval size, a 2-D case shows area, and a 3-D case shows volume. 


[0.9995,1] with S = 5e — 5. We can see in Table[I]that this ablation group has a 
propensity to reverse. The astute reader would notice that this trend does not 
seem to match Fig. [2j unlike the rest of our results. We have run simulations 
with the equations from |16j . and the simulations also produce reversal , not 
acceleration. The results of the simulation and model checking consistently dis¬ 
agree with the behavior denoted in Fig. [2] for this ALM group. We are currently 
investigating why this is the case. 


ALM, DVA Ablation Group: All three behaviors of this group are produced 
by varying only v 9 avm an d P 9 plm ■ Here reversal is satisfied over [0.02,0.0205] 


















with 6 = 2e — 5, acceleration over [0.67, 0.6705] with <5 = 5e — 5, and no response 
over [0.9995,1] with 5 = 5e—5. Repeating our experiments for this group, Table]!] 
shows the dominant reversal behavior. 

5.3 3-D Parameter Space 

Since the ablation groups we have used in this paper all feature at least one of 
the primary sensory neurons ( ALM , AVM , and PLM) ablated, we can only show 
the 3-D case for the original animal. 

For the 3-D case, in addition to p 9 ^y M and P 9 alm > we have the p 9 p[ M con¬ 
ductance. Finally, we get a non-zero value for no response in the control, but 
Table [l] shows that this value is an order of magnitude smaller than acceleration 
and several orders smaller than reversal. Reversal is satisfied over [0.01,0.0101] 
with S = 2e — 5, acceleration over [0.631, 0.6305] with <5 = 5e — 5, and no response 
over [0.63, 0.63005] with <5 = 5e — 5. 

6 Conclusions 

In this paper, we performed reachability analysis with discrepancy to automat¬ 
ically determine parameter ranges for three fundamental reactions by C. Ele- 
gans to tap-withdrawal stimulation: reversal of movement, acceleration, and no 
response. We followed the lead of the in vivo experimental results of m to ob¬ 
tain parameter-estimation results for gap-junction conductances for a number 
of neural-ablation groups. To the best of our knowledge, these results represent 
the first formal verification of a biologically realistic (nonlinear ODE) model of 
a neural circuit in a multicellular organism. 

Our results are further organized into how many of these three conductances 
we considered simultaneously. For each of these cases, we were able to determine 
parameter ranges for which the three TW responses are guaranteed to hold. 
Moreover, with the exception of the ALM- ablation group (an exception Wicks 
himself has noted about the ODE circuit model), our results match the relative- 
percentage trends of Fig. [2] 

Future work includes expanding the parameter ranges for TW responses, 
possibly by parallelizing the verification algorithm. We also plan to examine the 
additional ablation groups present in Fig. [2] 
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